function moments = generate_alt_cal_for_robustness(prefix_pf, prefix_skill, prefix, knitro_path)
    addpath('../')
    addpath('../Densities')
    addpath('../UtilityFunctions')
    addpath('../ProductionFunctions')
    addpath('../Economies/')
    addpath(knitro_path)

    dir_in = '../../../data/generated/matlab/Calibration/'
    dir_out = '../../../data/generated/matlab/simulation/'

    w0_range = 1e-4;
    wf_range = 1e4;

    % load calibrated parameters for skill distribution
    par_cal = tools.csv2struct([dir_in,prefix_skill,'_abilityDistThreeFixC_sol.csv']);
    par_mom = tools.csv2struct([dir_in,prefix_skill,'_abilityDistThreeFixC_mom.csv']);
    inequ_par = tools.csv2struct([dir_in,'inequ_par.csv']);
    constants = tools.csv2struct('../../../data/raw/constants.csv');

    DensSkill = abilityDistFixed(par_cal);

    DensPart = ParticipationNormal(par_cal);

    % utility function object
    par_util.epsilon = constants.epsilon;
    Util = UtilityQuasiLinear(par_util);

    % combine in SkillAndParticipation object
    par_dist.dens_skill = DensSkill;
    par_dist.dens_part = DensPart;
    par_dist.util = Util;

    par_dist.num_nodes_phi = 50;
    par_dist.scale_inc = par_cal.scale_inc; 

    par_dist.hsv_tax = true;
    Dist = SkillAndParticipation(par_dist);

    % Production function
    par_pf =  tools.csv2struct([dir_in,prefix_pf,'_KrusellRobotsEach_sol.csv'])
    mom_pf =  tools.csv2struct([dir_in,prefix_pf,'_KrusellRobotsEach_mom.csv'])

    %par_pf.eta_R = 0.9;
    
    par_pf.parse_inputs = true;
    par_pf.KrusellRobotsEach_use_properties = true;
    ProdFun = KrusellRobotsEach(par_pf);

    % Economy object
    par_econ.pf = ProdFun;
    par_econ.util = Util;

    par_econ.w_lb = 1e-6;
    par_econ.w_ub = 1e4;

    % robot prices and taxes
    inp.tau_B = par_pf.tau_B;
    inp.tau_E = par_pf.tau_E;
    inp.tau_S = par_pf.tau_S;

    par_econ.q_B = mom_pf.q_B;
    par_econ.q_E = mom_pf.q_E;
    par_econ.q_S = mom_pf.q_S;

    q_B_0 = mom_pf.q_B;

    par_econ.dist = Dist;

    par_econ.solver = 'knitro_eqsolve'

    par_econ.phi_lb = constants.phi_lb;
    par_econ.r = inequ_par.inequ_par;
    par_econ.par_dist = par_cal;

    par_econ.varnames = {'L_M','L_R','L_C', ...
                        'K_B_M','K_B_R','K_B_C', 'K_E','K_S',...
                        'tau_B','tau_E','tau_S','transfer','tau_hsv','lambda_hsv','compensation'};

    econ = KrusellWithRobotsEachEconomy(par_econ);
    econ.set_num_nodes(100,50,20,50);

    inp.L_M = par_pf.L_M_0;
    inp.L_R = par_pf.L_R_0;
    inp.L_C = par_pf.L_C_0;
    inp.K_B_M = par_pf.K_B_M_0;
    inp.K_B_R = par_pf.K_B_R_0;
    inp.K_B_C = par_pf.K_B_C_0;
    inp.K_E = par_pf.K_E_0;
    inp.K_S = par_pf.K_S_0;
    inp.transfer = par_cal.transfer;
    inp.tau_hsv = par_cal.tau_hsv;
    inp.lambda_hsv = par_cal.lambda_hsv;
    inp.compensation = 0;

    lb.L_M = 1e-16;
    lb.L_R = 1e-16;
    lb.L_C = 1e-16;
    lb.K_B_M = 1e-16;
    lb.K_B_R = 1e-16;
    lb.K_B_C = 1e-16;
    lb.K_E = 1e-16;
    lb.K_S = 1e-16;
    lb.tau_E = inp.tau_E;
    lb.tau_S = inp.tau_S;
    lb.transfer = inp.transfer;
    lb.tau_hsv = inp.tau_hsv;
    lb.lambda_hsv = inp.lambda_hsv;
    lb.compensation = 0;


    ub.L_M = 1e16;
    ub.L_R = 1e16;
    ub.L_C = 1e16;
    ub.K_B_M = 1e16;
    ub.K_B_R = 1e16;
    ub.K_B_C = 1e16;
    ub.K_E = 1e16;
    ub.K_S = 1e16;
    ub.tau_E = inp.tau_E;
    ub.tau_S = inp.tau_S;
    ub.transfer = inp.transfer;
    ub.tau_hsv = inp.tau_hsv;
    ub.lambda_hsv = inp.lambda_hsv;
    ub.compensation = 0;


    %% vesion with fixed capital (as in the calibration)
    econ.fix_other_capital = true;
    
    % also fix tax on robots
    inp.tau_B = par_pf.tau_B;
    lb.tau_B = inp.tau_B;
    ub.tau_B = inp.tau_B;

    % we want to compute the new equilibrium at the initial price
    droprange = [0];

    
    % compute baseline moments (corresponds to current calibration)
    moments = simulate(econ, inp, lb, ub, q_B_0, droprange, prefix, 'simulation_fix_capital', dir_out, par_cal, par_pf);        


function [results, results_percentiles] = simulate(econ,inp,lb,ub, q_B_0,droprange, prefix, suffix, dir_out, par, par_pf)
    sol = inp;
    stats = struct();
    percentiles = struct();
    results = table();
    results_percentiles = table();

    for i = 1:length(droprange)
        % change price of robots and recompute
        drop = droprange(i);
        econ.q_B = q_B_0 * (1-drop);
        
        fprintf('\n=================================================\n')
        fprintf('%s, price is %f\n', suffix, econ.q_B)
        fprintf('=================================================\n')
        
        if econ.fix_other_capital
            lb.K_E = inp.K_E;
            lb.K_S = inp.K_S;
            
            ub.K_E = inp.K_E;
            ub.K_S = inp.K_S;
        end
        
        tic
        [sol,fval,exitflag,output] = econ.compute_equilibrium(sol,lb, ...
                                                          ub);
        toc
        
        [stats, percentiles] = compute_stats(econ, sol, par, par_pf);
        stats.pricedrop = drop;
        stats.exitflag = exitflag;
        stats.fval = fval;
        results = [results; struct2table(stats)];
        results_percentiles = [results_percentiles; struct2table(percentiles)];
    end

end

function [stats, percentiles_strct] = compute_stats(econ, sol, par, par_pf)
    stats = sol;
    % assign prod fun pars
    stats.alpha = par_pf.alpha;
    stats.rho = par_pf.rho;
    stats.sigma = par_pf.sigma;
    stats.nu = par_pf.scale_L;
    stats.kappa_M = par_pf.kappa_M;
    stats.kappa_R = par_pf.kappa_R;
    stats.kappa_C = par_pf.kappa_C;
    stats.kappa_MR = par_pf.nu;
    stats.gamma_E = par_pf.lambda;
    stats.gamma_C = 1 - par_pf.mu;
    stats.gamma_B_M = 1 - par_pf.eta_M;
    stats.gamma_B_R = 1 - par_pf.eta_R;
    stats.gamma_B_C = 1 - par_pf.eta_C;
    stats.gamma_M = par_pf.upsilon;
    stats.A = par_pf.A;
    
    % assign other results
    stats.K_B = sol.K_B_M + sol.K_B_R + sol.K_B_C;
    stats.q_B = econ.q_B;
    stats.q_E = econ.q_E;
    stats.q_S = econ.q_S;

    [stats.Y_M,stats.Y_R,stats.Y_C,stats.Y_B_M,stats.Y_B_R,stats.Y_B_C,stats.Y_E,stats.Y_S]=econ.factor_prices(sol,par);
    stats.Y = econ.pf.Y(sol,par);

    w_lb = 1e-6;
    w_ub = 1e4;

    stats.part = econ.dist.mass_participants(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,par);
    stats.emp_M = econ.dist.emp_M(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.emp_R = econ.dist.emp_R(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.emp_C = econ.dist.emp_C(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);

    stats.avg_wage_M = econ.dist.average_wage_M(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.avg_wage_R = econ.dist.average_wage_R(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.avg_wage_C = econ.dist.average_wage_C(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);

    % welfare
    derived_vars = econ.compute_vars_from_input(sol);
    welfare = econ.welfare(derived_vars);
    stats.welfare_part_M = econ.welfare_participants_M(derived_vars);
    stats.welfare_part_R = econ.welfare_participants_R(derived_vars);
    stats.welfare_part_C = econ.welfare_participants_C(derived_vars);

    stats.welfare_non_part = econ.welfare_non_participants(derived_vars);
    stats.welfare_part = econ.welfare_participants(derived_vars);
    stats.welfare = econ.welfare(derived_vars);

    % tax revenues
    stats.tax_rev_inc = derived_vars.tax_rev_inc;
    stats.exp_transfer = derived_vars.exp_transfer;
    stats.tax_rev_B = derived_vars.tax_rev_B;
    stats.tax_rev_E = derived_vars.tax_rev_E;
    stats.tax_rev_S = derived_vars.tax_rev_S;
    stats.tax_rev_total = derived_vars.tax_rev_total;
    stats.excess_revenue = derived_vars.excess_revenue;

    % percentiles
    handle_h = @(w) 1/stats.part .* econ.dist.h(w,stats.Y_M, ...
                                                stats.Y_R,stats.Y_C,par);
    pct = 0.1:0.1:0.9;
    [percentiles, nodes, weights] = ...
        tools.Integration.compute_percentiles(handle_h,pct,1,200,100);

    % participation rates
    stats.part_M = econ.part_M(sol,par);
    stats.part_R = econ.part_R(sol,par);
    stats.part_C = econ.part_C(sol,par);

    % population shares
    integrand = @(w) econ.dist.dens_skill.f_M(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
    stats.pop_M = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);

    integrand = @(w) econ.dist.dens_skill.f_R(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
    stats.pop_R = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);

    integrand = @(w) econ.dist.dens_skill.f_C(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
    stats.pop_C = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);

    percentiles_strct.percentile = pct';
    percentiles_strct.value = percentiles;
    percentiles_strct.q_B = econ.q_B * ones(size(pct'));
    percentiles_strct.q_E = econ.q_E * ones(size(pct'));
    percentiles_strct.q_S = econ.q_S * ones(size(pct'));
    
    % income shares
    stats.labor_income_share = econ.labor_income_share(sol,par);
    stats.labor_income = econ.labor_income(sol,par);
    stats.M_share_in_labor_inc = (sol.L_M * stats.Y_M)/stats.labor_income;
    stats.R_share_in_labor_inc = (sol.L_R * stats.Y_R)/stats.labor_income;
    stats.robot_share_capital = (stats.q_B * stats.K_B)/(stats.q_B ...
                                                      * stats.K_B + ...
                                                      stats.q_E * stats.K_E + stats.q_S * stats.K_S);

end


end